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Abstract 

The self-modulation of waves propagating in nonlinear magnetic metamaterials is investigated. 
Considering the propagation of a modulated amplitude magnetic field in such a medium, we show 
that the self-modulation of the carrier wave leads to a spontaneous energy localization via the 
generation of localized envelope structures (envelope solitons), whose form and properties are 
discussed. These results are also supported by numerical calculations. 
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Increased research efforts have recently focused on metamaterials i.e., artificially struc- 
tured materials that exhibit electromagnetic (EM) properties not available in naturally oc- 
curing materials. For example, arrays of subwavelength split-ring resonators (SRRs) have 
been utilized, as proposed by Pendry et al. [1], to create magnetic metamaterials (MMs) 
with negative magnetic permeability /i up to Terahertz (THz) frequencies 

a a. 

There 

are only few natural materials responding magnetically at these frequencies. However, the 
magnetic effects in these materials are typically weak and often exhibit narrow bands (see 
Ref. Q and refs. therein). Thus, the realization of MMs at such frequencies will affect 
THz optics substantially, while it promises new device applications. Moreover, MMs with 
negative fi can be combined with plasmonic wires that exhibit a negative permittivity e J4], 
producing metamaterials which possess simultaneously negative values of fi and e 5]. The 
behavior of these metamaterials obeys a negative value of the refraction index, so they are 
usually referred to as negative index media or left-handed materials (LHMs). 

A concise analytical framework for the nonlinear behavior of LHMs was only quite re- 
cently proposed p, by considering the SRRs embedded in a Kerr-type medium, leading to a 
dependence of the effective values of e and fi on the EM field amplitude(s) of a propagating 
wave, roughly analogous to the Kerr effect in nonlinear optics. The SRR is modeled as a 
nonlinear RLC circuit, featuring an Ohmic resistance R, a self-inductance L, and a nonlinear 
capacitance C (due to the nonlinear dielectric filling its slit). Relying on this formulation, 
a nonlinear model for propagating modulated EM waves in LHMs predicts the possibility 
for modulational instability of EM waves and the occurrence of envelope bi-solitons, in fact 
exact solutions of a set of coupled nonlinear Schrodinger equations Q. The role of the 
anomalous self-steepening effect in modulational instability in LHMs has been studied by 
Wen et al. 8|. The spontaneous generation of a magnetic field due to the ponderomotive 
interaction of high-frequency EM fields with low-frequency electrostatic background oscilla- 
tions in the medium (plasmons) was suggested in Ref. [9]. A one-dimensional (ID) discrete 
array of nonlinear SRRs was considered in Ref. [loj |. which supports highly localized exci- 
tations in the form of discrete breathers. This model incorporates the nonlinearity of SRRs 
and the coupling among them at neighboring lattice sites via their mutual inductance M, 
thus accounting for Fourier dispersion of propagating charge waves. The combination of a 
nonlinear and dispersive lattice behavior allows one to anticipate the formation of nonlinear 
localized structures (solitons), eventually sustained by a balance among these mechanisms. 
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In this Brief Report, we consider the propagation of a weakly nonlinear charge variation in 
a nonlinear MM. Considering the self-interaction of the charge carrier wave, we show that 
the modulated wavepacket may be unstable, and that it may evolve towards the formation 
of localized structures in the form of envelope solitons. 

We consider an EM wave propagating in a ID array of N identical SRRs oriented in 
a direction perpendicular to the principal lattice axis, forming thus a magnetic-crystal- 
like arranges with lattice spacing D. I, the following we adopt the description (and 
notation) in Ref . [lCj , thus summarizing the essential building blocks of the theory in a self- 
contained manner, yet omitting unnecessary details. Considering a Kerr-type dependence 
of the permittivity e of the dielectric filling the SRR slits on the wave electric field E, viz. 
e = e (e; + a|E| 2 /.E 2 ), where e and ei respectively denote the permittivity of the vacuum and 
the linear part of the permittivity (E c is a characteristic electric field constant; a = +1/ — 1 
accounts for a focusing/ defocusing nonlinearity) , the electric charge Q n stored in the nth 
capacitor is given by Q n = CiU n [l + aU%/ (3e/?7 2 )], where Ci is the linear capacitance of the 
SRRs, U n = dgEg tn is the voltage across the slit of the nth SRR (viz. E 9 is the electric field 
induced along the slit) and U c = d g E c . 

Combining the above considerations, and approximating the nonlinearity by the first 

or expansion of its precise form (i.e., keeping up to cubic terms 



101]). the charge stored in the nth capacitor is found to obey the 



nonlinear term in a Tay 
in equation (6) in Ref. 
evolution equation 

-^{\q n -i -q n + Mn+i) ~q n + —q n = - s(t) , (1) 

where q n = Q n /(Cid g E c ) denotes the reduced charge, t = t/u>i = t(LCi) 1 ^ 2 , and A = M/L 
is the inter-site coupling constant. In the right-hand side (rhs), 7 = RCui accounts for 
the Ohmic impedance (friction) force F imp = jI n /(U c uJiCi) (J n = dQ n /dt) and e = £/U c 
is related to the electromotive force £ induced in each SRR due to the applied EM field. 
The rhs will be omitted in the following, i.e., by setting 7 = 5 = 0, thus namely neglecting 
Ohmic losses and electromotive forcing; these effects may be included a posteriori (e.g., as 
perturbations to the analytical solutions anticipated here) in future work. In what follows, k 
and uj represent their corresponding normalized values k — ► kD, and u — > u/ui. Then, the 
linear dispersion relation (DR) can be written in the form u = [1 — 2A cos(/c)]~ 1//2 , governing 



the propagation of magnetoinductive waves in such systems 



111 ] . This DR possesses a finite 
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cutoff at Umax = <u(k = 0) = 1 / \J 1 — 2A, and describes an inverse optic law, i.e., the 
group velocity v g = u'(k) = — Xu 3 sink is negative for all k's within the first Brillouin zone, 
viz. k e [0, 7t]. The wavepacket envelope therefore propagates (at the group velocity v g ) 
in the opposite direction with respect to the carrier wave (propagating at the phase speed 
v p h = u/k). The frequency band is therefore bounded by u min = u(k = 7r) = (1 + 2A) -1 / 2 
and uj max . Note that u max is increased by increasing the value of A, while u min is decreased 
in the same time. These results are apparently valid for A < 1/2, the physically meaningfull 
range for the system in consideration. A nonlinear generalization of the DR is obtained 
via a rotating wave approximation, i.e., by substituting q n = qexp[i(kn — ut)] + c.c. in the 
nonlinear form of Eq. ([1]) and retaining only first order harmonics. One thus obtains 

u 2 (k; \q\ 2 ) = (l - a\q\ 2 /e) (1 - 2A cos^))" 1 , (2) 

which incorporates the amplitude-dependence of the wave frequency. Assuming this depen- 
dence to be weak, and considering a modulated wave frequency u and wavenumber k close 
to the carrier values Uq and ko, respectively, one may expand as 



du 

u-u « - 



(k - ko) + 



2dk 

UJ=UJQ 



(k - k ) 2 + 



LU=Ll>Q 



d\q\ 2 



m 2 -\%\ 2 ), (3) 



9=90 



where go is a reference (harmonic wave, constant) amplitude. Considering slow space and 
time variables X and T, and thus setting uj — u —>■ id/dT and k — k —>■ —id/dX, one 
readily obtains the nonlinear Schrodinger- (NLS-)type equation 

where we have set tj) = q and ipo — Qo, an d defined the dispersion coefficient P = u"(k)/2, 
so that 

P = -Acj 5 (Acos 2 A; + cosA;-3A)/2, (5) 

and the nonlinearity coefficient Q = — (du / d\^\ 2 \^ = ^ o , viz. 

u a 1 u a 

Q= 2q1-||^o| 2 ^ 27r (6) 

Note that we have assumed \ipo\ 1, and thus neglected the dependence on \ip \ everywhere. 
Upon a Galilean transformation, viz. {X, T} —>■ {X — v g T,T} = {C) r } ; and a phase shift 
ip — > ipe~ l ^°^ T , one obtains the usual form of the NLS equation 

i£+i>g + WlV = 0, (7) 



which is known to occur in a variety of physical contexts [Id, Q- Eq- ©, along with 
the expressions (jSJ) and (jSJ), are the strong result of this calculation, to be retained in the 
analysis which follows. 

The evolution of a modulated wave whose amplitude is described by Eq. (j7j) essentially 
depends on the sign of the coefficients P and Q [13J. In specific, if PQ < the wavepacket 
is modulationally stable (and may propagate in the form of a dark- type envelope, i.e., a 
hole envelope soliton), while for PQ > the wavepacket is modulationally unstable. In the 
latter case, the wave may respond to random external perturbations (noise) by breaking- 
up to a "sea" of erratic oscillations (collapse) or (as suggested by analytical and numerical 
investigations 12, [l3, [ijj]) by localizing its energy via the formation of a series of localized 



envelope structures, i.e., propagating wavepackets modulated by a pulse-shaped envelope. 
We find that P is negative for low k, while it changes sign at some critical value k cr = 



cos -1 [(— 1 + \/T+ 12A 2 )/2A], thus acquiring positive values for k > k cr . On the other hand, 
the sign of Q is simply determined by the nature of the nonlinearity, i.e., Q is positive 
(negative) for a = +1 (-1). We may therefore distinguish two cases. For a = +1, the wave 
is modulationally stable (PQ < 0) for k < k cr , while it is unstable (PQ > 0) for k > k cr . 
For a = —1, the wave is modulationally unstable (PQ > 0) for k < k cr , while it is stable 
(PQ < 0) for k > k cr . 

To see this, first check that Eq. (J7|) supports the plane wave solution ip = 
ipo exp(iQ\ip \ 2 T); the standard linear analysis consists in perturbing the amplitude by set- 
ting: ^ = + 5 ^1,0 cos (k( — Cjt) (the perturbation wavenumber k and the frequency Co 
are distinguished from the carrier wave quantities, k and to). One thus obtains the perturba- 
tion DR Co 1 = Pk 2 (Pk 2 - 2Q|V>i,o| 2 )- O ne immediately sees that if PQ < 0, the amplitude 
ip will be stable to external perturbations. On the other hand, if PQ > 0, the amplitude 



ijj is unstable for k < J2Q/P\ipx j0 \; i.e. for perturbation wavelengths larger than a criti- 
cal value. The maximum perturbation growth rate is then found to be o~ max = \Q\ iV^ol 2 , 
and will therefore be inversely proportional to both e\ and E 2 [as can be seen by recov- 
ering dimensions in Eq. (jj^)]. This modulational instability mechanism is tantamount to 
the Benjamin- Feir instability in hydrodynamics, also long-known as an energy localization 
mechanism in solid state physics and nonlinear optics, among other physical contexts jl3 |. 
This type of analysis allows for a numerical investigation of the stability profile in terms of 
intrinsic medium properties. 
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Eq. (J7|) is integrable and possesses a number of exact solutions. Of particular interest to 
us are its constant profile, localized envelope solutions of the bright- (dark-) type, obtained 
for a positive (negative) value of PQ, which are of the form ip = ip Q exp(zO). The bright-type 



soliton solutions are given by 



15| 



e 



V'osech 



V L 



[v e (+(Q-v 2 e /2)r}/2P, 



(8) 
(9) 



where v e is the envelope velocity; L and Q represent the pulse's spatial width and oscillation 
frequency (at rest), respectively, and Ltp' = (2P/Q) 1//2 . Note that ip' is independent of the 
pulse (envelope) velocity v e here. 

A dark (black or grey) envelope wavepacket, looks like a propagating localized hole (a 
void) amidst a uniform wave energy region. The exact expression for black envelopes reads 
15]: 



e 



tanh(^ L 



— V P T 



[v e C + (2PQ^ 2 -v 2 e /2)r}/2P, 



(10) 
(11) 



where ip' L = (2\P /Q\) l l 2 . The grey-type envelope is given by 15] 



d sech' 



1/2 



(12) 



6 = [V C - (V 2 /2 - 2PQi)'\)T + e ]/2P - S sin 



-i 



</ tanh( — 



1 - rf 2 sech 2 ( i 



-1/2' 



Here, 6o is a constant phase, £ = £ — t> e r, and S 1 = sign(P) x sign(t> e — Vb). The pulse 
width L = \P / Q\ 1 / 2 / dipQ now also depends on the real parameter d, given by: d 2 = l + (v e — 



Vb) 2 / (2PQip'l) < 1. The (real) velocity parameter Vb satisfies the relation Vq — J2\PQ\ ip'l < 



v e < Vb + y2|PQ| ^'q . For c? = 1 one recovers the black envelope soliton. 

We performed numerical simulations using Eq. with the complete form of the nonlin- 



earity 



10j and initial conditions q n = tpo cos(n k — out), with i/jq in the form given in Eq. (jSJ) 



and Eq. fllOl) for bright and dark envelope solitons, respectively. The results are illustrated 
in Figs. 2 and 3 for both focusing and defocusing nonlinearity (a = +1 and —1), along with 
the corresponding analytic expressions, Eqs. ([HD and (fTOl . In Fig. 3, we have also depicted 
the envelope of the grey-type soliton, given by Eq. (Tl2|) with large "greyness", i.e., with 
d = 0.995. The numerically obtained envelopes are shown at some instant after they have 
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performed at least ten revolutions around the lattice, that is, more than 4200 time units 
(t.u.). However, they seem to be stable for much longer time intervals (at least up to 5 x 10 4 
t.u. for the cases we checked). This looks surprising, given that for the value of A used in 
the calculations (0.2 to be close to realistic values appearing in real systems) the system is 
still far from being continuous. This fact, along with the approximation of the nonlinear 
term, makes the expressions for the envelopes presented above to be the only approximate 
solution. Thus, they induce radiation in the lattice which eventually deforms the envelopes. 
Apparently, the simulations based on the discrete model show fairly good agreement with 
the analytical expressions derived from the continuous approximation of that model. In any 
case, the analytical expressions seem to reproduce, at least qualitatively, the observed en- 
velopes. We should note that our numerical simulation also provides evidence for dynamical 
dark soliton formation via modulational instability of a slightly modulated plane wave. This 
is a strongly nonlinear stage of the wave amplitude's evolution, which is not predicted by 
the linear amplitude perturbation theory employed above. The same procedure for bright 
solitons results in erratic oscillations. 

In conclusion, we have shown that the propagation of a modulated EM wave packet in a 
nonlinear MM, in the form of a lattice of SRRs, is characterized by amplitude modulation due 
to carrier wave self-interaction. An EM wave packet may be modulationally stable, and then 
propagate in the form of a localized void (a hole, amidst a uniform charge density), or it may 
be intrinsically unstable, and thus possibly evolve towards the formation of envelope pulses 
(bright solitons). Explicit expressions for these nonlinear magnetoinductive excitations, 
whose existence is supported by numerical calculations, are provided in terms of the intrinsic 
material parameters. These results are of relevance in metamaterial-related applications, in 
materials which may be "tuned" appropriately in order for the forementioned excitations to 
occur, pretty much like optical fiber solitons in nonlinear optics. 
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Figure captions 



Figure 1 

(a) The linear dispersion u = oo(k), and (b) the dispersion coefficient P = P(k), for 
A = 0.05 (black-solid curve), 0.2 (red-dotted curve), 0.35 (green-dashed curve). 

Figure 2 

Dark-type modulated wavepackets, for A = 0.20, ee = 2, ip' Q = 0.1, iV = 100, and (a) 
a — —1, k — 1.12 > k cr (u = 1.10), L = 3.74; (b) a — +1, k — 0.86 < k cr (u = 1.16), 
L = 4.51. The filled circles correspond to the g n 's, with q n the charge at site number n, while 
the black-dotted curve serves as a guide to the eye. The red-solid curves are the envelope 
Eq. (flUj) . and the green-dashed curves the envelope Eq. (fl2l) . with parameters as in (a) and 
(b), and d = 0.995. 

Figure 3 

Bright-type modulated wavepackets, for A = 0.20, ee = 2, ip' = 0.1, iV = 100, and (a) 
a — —1, k — 0.86 < k cr (u = 1.16), L = 9.02; (b) a — +1, k — 1.12 > k cr (u = 1.10), 
L = 7.48. The filled circles correspond to the g n 's, with q n the charge at site number n, while 
the black-dotted curves serve as a guide to the eye. The red-solid curves are the envelope 
Eq. ([8]), with parameters as in (a) and (b). 
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FIG. 1: (color online), (a) The linear dispersion to = u>(k), and (b) the dispersion coefficient 
P = P(k), for A = 0.05 (black-solid curve), 0.2 (red-dotted curve), 0.35 (green-dashed curve). 
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FIG. 2: (color online). Dark-type modulated wavepackets, for A = 0.20, ei = 2, tfj' = 0.1, N = 100, 
and (a) a = —1, k = 1.12 > k cr (uj = 1.10), L = 3.74; (b) a = +1, k = 0.86 < k cr (to = 1.16), 
L = 4.51. The filled circles correspond to the g n 's, with q n the charge at site number n, while the 
black-dotted curve serves as a guide to the eye. The red-solid curves are the envelope Eq. (fTUj). and 
the green-dashed curves the envelope Eq. (fT2"j) . with parameters as in (a) and (b), and d = 0.995. 
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FIG. 3: (color online). Bright-type modulated wavepackets, for A = 0.20, = 2, ip' = 0.1, 
N = 100, and (a) a = -1, = 0.86 < /c cr (cj = 1.16), L = 9.02; (b) a = +1, k = 1.12 > fc cr 
(a; = 1.10), L = 7.48. The filled circles correspond to the g n 's, with q n the charge at site number 
n, while the black-dotted curves serve as a guide to the eye. The red-solid curves are the envelope 
Eq. (jSJ), with parameters as in (a) and (b). 
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